quadmesh

Build a “quadmesh”" in R.

library(quadmesh)
library(raster)

## Loading required package: sp

data(volcano)
r <- setExtent(raster(volcano), extent(0, 100, 0, 200))


qm <- quadmesh(r)

library(rgl)

## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display

## Warning: 'rgl_init' failed, running with rgl.useNULL = TRUE

scl <- function(x) (x - min(x))/diff(range(x))
shade3d(qm, col = grey(scl(qm$vb[3,qm$ib])))

rglwidget()

A “quadmesh” is a dense mesh describing a topologically continuous surface of 4-corner primitives. I.e. a grid, without the “regular”. This is useful particularly when combined with map projections and texture mapping.

We are not limited to a regular grid, trivially let’s distort the mesh by a weird scaling factor.

The topology of the grid is still sound, but we are no longer bound to the regular constraint.

qm <- quadmesh(r)

qm$vb[1,] <- qm$vb[1,] * qm$vb[2,]/54
open3d()

## null 
##    3

shade3d(qm, col = grey(scl(qm$vb[3,qm$ib])))

rglwidget()

Why quadmesh?

Why meshes at all?

The simplest kind of mesh is a basic raster. Consider a matrix of 12 values.

(m <- matrix(1:12, nrow = 3))

##      [,1] [,2] [,3] [,4]
## [1,]    1    4    7   10
## [2,]    2    5    8   11
## [3,]    3    6    9   12

On its own this matrix has absolutely nothing to do with spatial data, it is literally a collection of 12 numeric values in a given order, and by the magic of programming we’ve nominated a shape of 3x4. We can’t help but think about this shape spatially however, but there’s a problem. Does each element occupy space or should we consider them to be infinitesimal locations?

R provides either interpretation (to simplify this story we nominate locations for the rows and columns explicitly).

When considered as an image, each matrix element occupies a certain space in width and height, but when considered as a point set the numbers simply float at the given locations. Which is correct? (Spoiler: Both are correct, it simply depends what we are doing.)

x <- seq(1, nrow(m)) - 0.5
y <- seq(1, ncol(m)) - 0.5
image(x, y, m)

text(expand.grid(x, y), lab = m[])

The raster package defaults to the image interpretation and helpfully assumes the values are nominally at the centre points as shown above. We have to nominate the extent or we end up in 0,1 range, we also have to invert the order of the values because raster counts from the top of the page and R’s matrix uses column-major order.

library(raster)
(r <- raster(t(m[, ncol(m):1]), xmn = 0, xmx =ncol(m), ymn = 0, ymx = nrow(m)))

## class       : RasterLayer 
## dimensions  : 4, 3, 12  (nrow, ncol, ncell)
## resolution  : 1.333333, 0.75  (x, y)
## extent      : 0, 4, 0, 3  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 1, 12  (min, max)

R’s image and rasters in general are so efficient because they only store this minimal amount of information: the actual data values, and the extent and dimensions of the space they occur in. If we had to store the centre coordinate of every cell, or worse the corner coordinates then the data storage goes up dramatically. Every software that deals well with these kinds of data has to treat these coordinates as implicit. We can easily expand the centre coordinates.

xyz <- as.data.frame(r, xy = TRUE)
head(xyz)

##           x     y layer
## 1 0.6666667 2.625    10
## 2 2.0000000 2.625    11
## 3 3.3333333 2.625    12
## 4 0.6666667 1.875     7
## 5 2.0000000 1.875     8
## 6 3.3333333 1.875     9

tail(xyz)

##            x     y layer
## 7  0.6666667 1.125     4
## 8  2.0000000 1.125     5
## 9  3.3333333 1.125     6
## 10 0.6666667 0.375     1
## 11 2.0000000 0.375     2
## 12 3.3333333 0.375     3

but to expand the corners we have to jump through some hoops and even then we get every instance of the corners, not only for each cell but to explicitly close the cell as a polygon.

as(as(raster::rasterToPolygons(r), "SpatialLinesDataFrame"), 
   "SpatialPointsDataFrame")

## class       : SpatialPointsDataFrame 
## features    : 60 
## extent      : 0, 4, 0, 3  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 4
## names       : layer, Lines.NR, Lines.ID, Line.NR 
## min values  :     1,        1,        1,       1 
## max values  :    12,       12,        9,       1

There are only 20 unique coordinates at the corners, which is where quadmesh comes in.

library(quadmesh)
qm <- quadmesh(r)
str(qm)

## List of 6
##  $ vb           : num [1:4, 1:20] 1.49e-08 3.00 1.10e+01 1.00 1.33 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:4] "x" "y" "z" "1"
##   .. ..$ : NULL
##  $ ib           : int [1:4, 1:12] 1 2 6 5 2 3 7 6 3 4 ...
##  $ primitivetype: chr "quad"
##  $ material     : list()
##  $ normals      : NULL
##  $ texcoords    : NULL
##  - attr(*, "class")= chr [1:2] "mesh3d" "shape3d"